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THE GODDARD VERSION OF THE SCHUBART-STUMPFF N-BODY PROGRAM 


P. A. Cornelia 
B. E. Lowrey 


Abstract 


The S chub art -Stumpff N-Body Program computes solar system orbits 
of the planets and of bodies of zero mass. It can also be used to 
solve more general problems in mechanics. A Fortran IV adaptation of 
the Yale version of the program is available to users from the GSFC 
documentation center, Primary modifications ' to the Yale program in- 
clude simplification of program input and the addition of an option 
to compute the osculating orbital elements. The present document 
briefly summarizes the major features of the program and the 
development of the Schubart-Stumpff initial values, discusses the 
new input requirements and modifications, and presents a sample case 
for clarification. 
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X. Method of Schubart -Stumpf f Program 

S chub art and Stumpf f (1966) have chosen to regard the orbital 
computation problem as one of finding the solution to a system of 
simultaneous differential equations, Basic to this approach is the 
decision to treat all bodies alike rather than to use the special 
perturbation approach which distinguishes perturbed from perturbing 
bodies. Hence only initial values and an appropriately chosen step- 
size are needed to achieve solution. The method of integration used 
is that of Adams - St ormer: a difference method with constant step-size 

for the numerical integration. 

By avoiding the special perturbation method it is unnecessary to 
input coordinates of perturbing bodies at other than the starting 
epoch, thus avoiding much of the data manipulation associated with 
reading in tables of perturbing bodies, which of itself can be a 
formidable problem, both logically and logistically. 

The method of numerically integrating the massive bodies rather 
than using tabular input is not uncompetitive in machine time, depending 
upon the problem. If it is desired to compute the orbits of many 
massless bodies for the same period Of time, the machine time required 
to compute the planetary orbits becomes a small fraction of the total 
time. But the N-Body Program also has the ability to study the orbits 
of the planets themselves, as for example, Lieske's (1967) preparation 
of JPL's Development Ephemeris Humber 28 . In another application, 
Marsden (1969) used the program (with some modifications for 


- 1 - 



differential corrections) to determine the influence of non-gravitational 
forces on the orbits of short period comets. 

It is important to note that there are no provisions made for the 
problems that arise in extremely close passages. The planets are taken 
as mass points. Also the experimenter must know at what epoch a close 
approach will occur since there is no automatic control of step-length. 
Then he may interrupt the computation, input a smaller step-size with 
which to compute for the duration of the approach, and then again 
interrupt after approach to input a larger step-size. Schubart and 
Stumpff chose this seemingly clumsy method of step-length control 
because in their own theoretical work they could predict the epoch of 
a close approach very easily. The authors of the present note have 
retained this method in the GSFC version of the program. They plan, 
however,, to implement a Nordsieck-type of predictor-corrector method 
in order to permit dynamic internal calculation of optimum interval- 
size for each integration step. 

II. The Initial Conditions 

The Schubart-Stumpff N-Body Program together with the proper 

t 

initial conditions for the planets calculates to a high degree of 
accuracy orbits of bodies of (essentially) zero mass and simultaneously 
the orbits of the planets. 

For solar system orbits Schubart and Stumpff have derived a set 
of starting conditions of the planets for epoch JD = 2430000.5 . These 
derivations are discussed in great detail in their paper. The 
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following paragraphs summarize their work. 

With their initial values and a step-size of 2 days, the program 
reproduces to 10 decimal places the ephemerides of the planets of the 
solar system under the following conditions: 

(1) Relativistic effects are ignored; 

(2) The perturbing effects of Mercury are only approximated: 
Mercury's mass is added to that of the sun, thus introducing an error 
of 10 7 A.U. in the location of the origin of a heliocentric system; 

(5) The mass of the moon is added to that of the earth hut the 
perturbations in the earth-moon orbit caused by the moon are not con- 
sidered. 

(4) Perturbations caused by Pluto are ignored, a decision 
motivated not only by economics but also by the fact that Pluto's mass 
is not known with sufficient accuracy. 

The Eckert, Brouwer, Clemence ephemeris (1951) of the five outer 
planets (Jupiter to Pluto); Herget's computations from Newcomb's 
tables (1953, 1955) of the ephemeris for Venus and the center of mass 
of the earth-moon system; and the R. L. Buncombe -G. M. Clemence 
ephemeris (i960, 1964) for Mars were the standards used in the de- 
velopment of the initial conditions and in the comparison of the final 
results. 

To facilitate conversions and comparisons 11th differences of the 
acceleration were used which correspond to the value of 5 for M in the 
Fortran code. The initial epoch, J.D. = 2430000.5 > was chosen because 
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for periods of 400 days on either side of that date, the Eckert, 

Brouwer, Clemence observations of the 5 outer planets were especially 
good, a factor vital to deriving the starting velocities. 

For any solar system problem the Schubart-Stumpff initial values 
can always be used. If the epoch for which the experimenter wishes 
to enter coordinates of zero mass bodies differs from their epoch, 
a negative step-size may be used to integrate back to an earlier epoch, 
a positive step-size to integrate forward to a later epoch. The 
calculations then proceed with the additional bodies. 

III. GSFC Version of the Schubart-Stumpff N-Body Program 
A. Modifications 

The original Schubart-Stumpff N-Body Program was written in the 
mid-1960's in the Fortran II language for use on a smallish computer 
at the University of Heidelberg. Subsequently, Schubart and M. Cooke 
( 1965 ) rewrote the program in the Fortran IV language for use on an 
IBM 709^ at Yale. It is the Yale version that the authors have modified 
for use on Goddard Space Flight Center's IBM S/ 36 O computers. 

These modifications are simplifications in the authors ' opinion: 
the Yale version still reflected the idiocyncrasies of the Fortran II 
language in coping with double word computations and complicated data 
input. These peculiarities made it clumsy to input data to the 
program. The authors adopted the NAMELIST convention of the Fortran IV 
language to input all data, except for one initialization card. 

The NAMELIST option, by allowing selective initialization of data 
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without destruction of data in locations not specifically named on a 
given READ, enables the programmer to eliminate much coding of trigger 
recognition parameters and statements. Thus although the input is 
still flexible, only one read statement is required for the entire 
program to initialize all (but two) of the input parameters. 

A second modification was the addition of an option to print the 
osculating elements. A parameter triggers a call to a subroutine 
written by Blanchard and Wolf ( 1967 ) which converts the position and 
velocity vectors at a given time into the orbital elements by means 
of the two-body formulas. 

Finally, the addition of 3 parameters permits the positional co- 
ordinates, the velocity coordinates and the mass parameters to be input 
with incompatible dimensional systems, by converting all to a compatible 
system following input. Hopefully this will eliminate the tedious 
hand computations which might be required for some problems. 

B. Subprograms of the Program 

1. MAIN - determines size of the A, B, and D [large] arrays, 
which initial conditions are to be used, and calls CONTRL. 

2. SUBROUTINE CONTRL - drives the program by processing in- 
put, calling integration routines and requesting output. 

3. SUBROUTINE KOEFFZ - calculates the K-coefficients, (CAl), 
used in the starting iteration. 

4. SUBROUTINE ANFITN - performs the starting iteration and 
converts, if necessary, positional and velocity beginning co-ordinates 
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from a system originating in body #1, to a barycentric system, used, 
in the integration. 

5. SUBROUTINE WW - computes the backward differences of the 
accelerations, (BESCHL) . 

6. SUBROUTINE SCHRTT - calculates from BESCHL, the backward 
differences of the co-ordinates (XNABLA) and the co-ordinates (x), 
thus integrating from time T to time T + DELTAT. 

7. SUBROUTINE DRUCKE - controls output at constant time 

intervals . 

8. SUBRQT JTIHE RT.EMNT - controls output at arbitrary time 

intervals * 

9. SUBROUTINE ORBIT - converts position and velocity co- 
ordinates (for output purposes only) into osculating elements at time 
intervals specified by DRUCKE or ELEMNT. 

C , Input Parameters 

Immediately following the IBM S/36O JCL card, // G0 . DATA5 DD* 
there is one card required. Then the NAMELIST data set(s) with para- 
meters initialized (in arbitrary order) follows. 

1. Card No. 1 : This first card sets values of 2 parameters 

NSIZE, occupying columns 1-5 of the card, right -adjusted, and INIT, 
occupying columns 6-10, right -adjusted. 

NSIZE determines the sizes of 3 large arrays — A,B,D — these sizes 
being dependent upon the number of bodies in a run of the program: 

50 s NSIZE a N. By setting this parameter himself, the programmer is 
able to reduce the memory requirements of the program, which may mean 
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a higher priority in an MVT environment, hence a faster turn-around 
time. 

HUT. The program initializes the coordinates and mass parameters 
for the Sun-Mercury system and the planets Venus to Pluto to those 
values that Schubart and Stumpff derived for J.D. = 2430000. 5 . It 
also sets the integration parameters to the values that Schubart and 
Stumpff considered optimum. 

By setting INIT = 0, the programmer can use these initials con- 
ditions. The NAMELIST data set is then used to initialize print -punch 
parameters, enter coordinates of massless bodies and alter any pre- 
initialized parameters (if he -wishes). 

By setting UTIT = 1, the Schubart-Stumpff conditions are over- 
ridden. Hence all pertinent initializations must be made in the 
NAMELIST data set . 

2. NAMELIST parameters : The first card of the NAMELIST 

DATA SET must contain & INPUT, with the & in column 2. Subsequent 
cards start in column 2; & END terminates the data set. Those users 
unfamiliar with the NAMELIST convention of the FORTRAN IV language 
are referred to the IBM manual on the FORTRAN IV language. 
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N number of bodies being integrated. As new bodies are added 

during a computation at an epoch different from the starting 
epoch, W must be increased accordingly. 

1 i NSIZE 

EM(l) mass of the ith body in units of m. If a relative system is 
used, the origin must be EM( l) : 

KQ = k 2 = G # EM( l) where G is the universal gravitational constant 

in units of L' 3 /S' 2 m, where L' is the unit of length, S' is the 
unit of time, EM(l) is the mass of body #1 in units of m. 

W is the conversion factor such that WW expresses KQ in units of 

L 3 /S 2 m, where L & S are the units of integration, 

XP(l,l) The x,y,z components of the position of ith body in units of, L' 1 , 

XP(2,I) 

XP(3,l) in the appropriate coordinate system (barycentric or relative to 
body #l) . 

DIST is the conversion factor necessary to express XP in units of 

L. DIST. = L/L" 

XDOT(l,l) the u,v,w components of the velocity of the ith body in units of 

XDOT( 2,l) 

XDOT( 3,l) L" ' /S' " , in the appropriate coordinate system. 

VEL is the conversion factor necessary to express XDOT in units of 

L/S . VEL = L/s/L ' " /S " ' . 

H is the integration step length in units of time S. 

T is the starting epoch for the problem in arbitrary time units . 

This is the time that appears in the output , 
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DELTAT is the integration step- length in the same units as T. If con- 
verted to units of S, DELTAT is equal to H. At the start of each 
integration step T is incremented by DELTAT. 

M order of the integration (lVK5) in the initial iteration. This 

corresponds to 2#M differences; in the extrapolation components 
to 2*M+1 . M=5 is the value used by Schubart and Stumpff in their • 

calculations . 

IEG=1 Input coordinate system has its origin in body # 1. 

IEG=0 Input coordinate system is barycentric, the coordinate system 

used in the integration. Whenever IEG = 1, the position and 
velocity components are converted from a relative to a barycentric 
system. 

IEXP 10*** ( -IEXP) is the limit of accuracy for the initial iteration. 

Note : The Schubart-Stumpff initial conditions for solar system orbital 

computations are coded into the program. These involve the initialization 
of the following parameters: 

EM(l), EM( 2) , ..., EM(9), 

XP( K, I) , , XP(K,9), K=l,2,3 

XD0T(K,l), , XD0T(K,9) K=l,2,3; 

the parameters for the Sun-Mercury System, Venus , Earth- Moon System, Ifers, 
Jupiter, Saturn, Uranus, Neptune and Pluto, respectively; as well as the 
parameters 

N, W, KQ, BIST, VEL, H, M, IEG, IEXP, T, DELTAT 
IORB activates calls to SUBROUTINE ORBIT, which computes osculating 

orbital elements, from subroutines CONTRL,ELEMNT,DRUCKE and 
prints them as specified by values of IORBIT. 
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I0RB=0 


I0RB=1 

iorbit(j) 

iorbit(j) 

N7 


F7=0 

N7=k 

ILEM(l) ,T- 


KD 


KD=1 

KD=2 

KD=3 


do not compute any orbital elements, i.e. do not call Orbit, 
compute orbital elements fear selected bodies. If I0RB=1, then 
it is necessary to specify IORBIT(j), j»l, . ..,N. 

•0 Do not compute or print osculating orbital elements for 
body # J. 

=1 Compute for body #J, and print the values. 

determines whether or not SUBROUTINE CONTROL initiates calls to 
SUBROUTINE ELEMNT, a print-punch control routine. SUBROUTINE 
ELEMNT is used when output at arbitrary integration intervals 
is desired. 

No calls to SUBROUTINE ELEMNT 
k calls to SUBROUTINE ELEMNT 

=1, . . . ,N7 Associated with each integration step is a step number, 
the ILEM array specifies in increasing order at which step 
numbers calls to SUBROUTINE ELEMNT are to be made. 

Note: ILEM(l)^M, else SUBROUTINE ELEMNT is never called even 

though N7>0. 

is a parameter which is utilized in SUBROUTINE ELEMNT to determine 
which coordinate system to use in the output and the number of sig- 
nificant figures to print. For each body in the problem the co- 
ordinates are printed. Error bound information is also supplied 
as well as time, T . 

Coordinate system relative to body #1. 

Single precision. 

Coordinate system: barycentric . 

Single precision. 

Coordinate system relative to body #1. 

Double precision. 
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KD=k 


IED 

IED=KD 

IED=KD+4 

IZA, IDELTZ 
IZA 

IDELTZ 

IZN 


ifunch(j) 


ifunch(j)= 


Barycentric system 
Double precision, 

A control parameter used in SUBROUTINE ELEMNT to specify type 
of output. 

Print output form using ELEMNT in form determined by value of KD. 
Punch or write-on tape, as well as print. The Job Control Language 
(JCL) Statements will direct whether to use the punch or tape. 

,IZN initiate calls from SUBROUTINE CONTRL to SUBROUTINE DRUCKE. 
first step number at which DRUCKE called. 

Every IDELTZ step after IZA SUBROUTINE DRUCKE is called until 
step IZN reached, 

SUBROUTINE DRUCKE called. This is the last iteration step as 
well: the program, depending upon the value of IG0, does one 
of the following: terminates the run, begins a new case, or 
continues the case by adding bodies of zero mass or changing the 
integration step-size. 

controls informtion written concerning jth body in SUBROUTINE 
DRUCKE. This information is as follows: step number, time, co- 

ordinates as well as error growth information. DRUCKE allows 
selection of bodies for which to print /punch whereas ELEMNT 
print s/punch for all bodies. 

1 Barycentric system. Double Precision Print. 

2 Barycentric system, D.P., Punch 

3 Barycentric system, D.P., Print and Punch 
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=4 Bary centric System, f P., Print 
=5 Barycentric, S.F., Punch 
=6 Barycentr: c, S.P., Print and Punch 
=7 Relative, D.P., Print 
=8 Relative D.P., Punch 
=9 Relative D.P. Print and Punch 
=10 Relative S.P., Print 
=11 Relative S.P., Punch 
=12 Relative S.P., Print and Punch 
=15 No output for Body #J. 

When IPUNCH( J)>14 it signals to the CONTRL Program that Body #J is of zero 
mass and has been added to the program at an epoch different from the 
starting epoch. 

The output options for Body #J then are such that 14 & IPUNCH(j) jS 26 and 
such that IPUNCH(j)-15 gives the correct option. 

IGO is a SUBROUTINE CONTROL parameter which tells the program what 

to do when Step IZN is reached ( IG0=1 or IG0=2) or what tests to 
make for incorrect data on a continuation case (lG0=5) • 

IG0=1 when IZ=IZN go to 2 and continue case by 

(1) adding new bodies of zero mass and incrementing N accordingly 

and/or 

(2) decreasing H and BELTAT for a close approach and modifying 
IZA,IZN and IDELTZ if necessary or 

( 3 ) increasing H and DELTAT when close approach calculations 
completed and modifying IZA,IZN and IDELTZ if necessary. 
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IG0=2 


when IZ=IZN, go to 100: Initialize the NAMELIST DATA SET and 

begin new case. 

IGO-3 is used in a continuation case, where if new bodies are added, 
they are checked to see that they are of zero mass with proper 
coordinates, relative to body #1. 

D. Internal Conversion of Input 

1. Units 

The position and/or velocity components may be converted internally 
to the units required for computation using the DIST and/or VEL parameters, 
respectively. W may be used to convert the force, KQ. 

If KQ=k 2 , where k is the Gaussian constant for solar system bodies, 
k®. 017202 098 950, then the units are distance in a.u., time in years, 
mass of sun = 1. Schubart and Stumpff have used a unit system where 
velocities are in a.u./ifO ephemeris days; hence they set KQ^O^ 2 = 

.474 3^5 981 216 687, EM(l) = 1. (The starting values for INIT=0 reflect 
this) . If dimensionless mass units are not desired, the force unit, KQ, 
may be adjusted instead. Whenever it is not necessary to convert distance 
units, set DIST=1; velocity units, set VEL=1; mass units, set W=l. 

2. Coordinate System 

If IEG=1, the program assumes that the origin of the relative 
coordinate system used to input the positional and velocity components 
is body #1. It then converts to a barycentric system which is used for 
computation. The position and velocity coordinates of body #1 must be 
set equal to zero if the relative system is used. 
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E. Program Notes 

1. Output Parameters 

Associated with each integration step, which is of length, H, are 
two parameters, T and IZ. T is the time in arbitrary units and is 
initialized in the NAMELIST DATA SET to the starting epoch. It is in- 
cremented by DELTAT, also initialized in the NAMELIST set . DELTAT, if 
converted to units of time, S, used in the integration, is equal to H. 

T is the parameter which is printed/punched in the output whenever time 
is specified. IZ is the step-number. It is initialized internally and is 
incremented by 1 at the start of each integration step. The ILEM array 
is compared to IZ and whenever an element of ILEM equals IZ, a call to 
SUBROUTINE ELEMNT is made for printed or punched output according to the 
values of KD and IED for all N bodies . This comparison enables output 
at arbitrary step numbers. The IZA,IBELTZ, and IZN parameters are used 
to get output at regular step intervals using SUBROUTINE DRUCKE. In 
this subroutine the IPUNCH array specifies the type of output desired for 
each individual body. 

Whenever I0RB=1, SUBROUTINE ELEMNT and SUBROUTINE DRUCKE call 
SUBROUTINE ORBIT to compute and print the osculating elements for bodies 
as specified in the IORBIT array. 

2. Machine Time: The running time is proportional to the number 

of distances which must computed at every integration step. There are 
nx(n! - l)/2 + ni*n 0 such distances where nj. is the number of bodies 
with mass and n 0 is the number of bodies of zero mass. Thus, when one 
computation has been obtained for a particular computer, the machine time 
can be predicted for other computations on that computer. 
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IV. Sample Case - Orbit of Lost City Meteorite 

The Lost City Meteorite struck the earth on January 3, 1970* The 
Prairie Network observed its entry into the earth’s atmosphere and 
provided the data which were used in the N-Body program to compute its 
(probable) heliocentric orbits. 

The computations were done in two parts. Using the Schubart-Stumpff 
initial conditions for Julian Day 2430000.5, the integration was first 
carried forward (H > 0 and DELTAT > 0) to Julian Day 2440590, 

(January 3, 1970 ), i n order to provide the planetary coordinates at 
the time the meteor was observed. These coordinates plus the observed 
coordinates of the meteor then provided the initial conditions for an 
integration backward (H < 0 and DELTAT < 0) through time for 300 years. 

The input data for the first part are shown in Figure 1. The 
first line sets the dimensions (20) of the A, B, and D matrices to 
A(20,3, 12), B(20,3,12), D(20,20), respectively, in the subroutines in 
which they appear corresponding to their initializations to the same 
sizes in the EEAL*8 statement of the MAIN program. 

The second number (l) instructs CONTRL to override the programmed 
initializations . 

The following discussion will use the left most numbers (line 
numbers) in referencing parameters. 

A. Line 100: 

1. IEG = 1: The Schubart-Stumpff planetary coordinates are 

relative to Body TL, ( Sun -Mercury system). 
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B . Line 200 : 


1. KD = 3: when in SUBROUTINE ELEMNT print coordinates 

in relative system in double precision. 

2. IED = 3=KD: suppresses punching of the coordinates when 

they are printed. 

C. Lines 300-2000: 

These lines contain the Schubart-Stumpff planetary coordinates for 


J.D. 2430000.5 with XP(1,K), XP(2,K), XP('3,K), K = 1,2,... 9 , 

the positional coordinates relative to body #1 of the body in 
a.u. and XD0T(l,K), XD0T(2 ,k), XD0T(3,K), K = 1,2, .... 9 , the 


velocity coordinates relative to body ^1 of the body in a.u./40 
ephemeris days. 


K = 1 

Sun-Mercury system 

K = 2 

Venus 

II 

VM 

Earth-Moon system 

K = 4 

Mars 

K = 5 

Jupiter 

K = 6 

Saturn 

w 

II 

— v] 

Uranus 

K = 8 

Neptune 

II 

VO 

Pluto 

Lines 

2100 - 2200 : 


Mass parameters, normalized to the sun of bodies, K = 1,2. ...... . 9 . 
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E. Lines 2300-2800: 

1. H = 9 the number of bodies integrated in this first part. 

2. KQ = .473^5961216687, reflecting the fact that the 
velocity units are a.u./40 days. 

3. W = 1 because the computation of KQ was done externally; 
had it been done internally, one could have used KQ =.017202098950 

and W = 40 to the same end. 

4. Because the input units for the coordinates and the 
integration units were compatible, the program initialized values of 
DIST and VEL did not have to be changed. 

5 . T is set arbitrarily to zero although is could have been 
initialized to a (more) meaningful value, e.g., T = 2430000.5 • 

6. Since the unit of time is 40 ephemeris days, the integration 
step-size, H, of 2 days is set to .05 (of 40 days) while the printing 
parameter, DELIAT, is correspondingly set to 2 days. 

7 . M = 5 determines the orders of the predictor -corrector 

method. 

8. N7 = 5 determines that CONTRL will CALL ELEMNT for out- 
put five times during the run at steps numbers given by the ILEM array: 
the lO^ 1 , 30 t5a , 500^, lOOO^ 1 , and 1059^ steps. ILEM(l) must always 
by greater than or equal to M . 

9 . IZA., IDELTZ, IZH signal COMTRL to call the output sub- 
routine, DRUCKE, starting at step number 1059, every 1059^ step 
thereafter, ending at step number 5295 , which is also the last 
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integration step. At this time IGO = 1 tells CONTRL to return to 
statement number 2 for additional information to continue the case. 

10. IQRB = 1 instructs ELEMNT and DRUCKE to CALL ORBIT at 
each output step with orbital elements being printed for bodies 2 

through 9, IORBIT(k) = 1, K = 2,3,4, ,9 but not for body #1, 

I0EBIT (l) =0 . 

11. IFUNCH =9*7, 41*13 instructs DRUCKE to print in double 
precision the planetary coordinates relative to body #1 of bodies 1 
through 9» 

F. Note that the namelist begins at line 100 with &HTPUT, 
starting in column 2 or greater and concludes at line 2800 with &END. 

G. The second part (Figure 2) of the run begins with the input 
of the position and velocity coordinates of the Lost City Meteorite 
[xp(j,10),j =1,3 a n<i XD0T(j,10),j = 1,3, respectively], and the cor- 
responding change of N to 10. To integrate backward, H becomes -0.5 
and DELTAT, -2 . For printing T is re-initialized to 2440590. In 
order to be consistent with the astronomical convention of standard 
days, IZA, IDELTZ, and IZN are re-set so that printing (in DRUCKE) 
occurs in multiples of 400 days, starting at T = 2430000.5 , with IZN 
such that the computations continue back through 300 years. IORBIT(l) = 1 
and IFUNCH(lO) = 7, enable output of the coordinates (DRUCKE) and 
orbital elements (ORBIT) of the Lost City Meteorite. IGO = 3 
indicates that the ensuing computations are a continuation of the run 

and not the start of a new case. 
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MAIN PROGRAM SETS HI MENS ION SIZES ON* A f B , A NO 0 MATRICES. 
NSIZE... ... SFTS MAXIMUM MIIMRFR OF BODIES FOR A RUM. 

INI T* . * • * • * INDICATES WHETHER USER WILL INPUT HIS OWN INITIAL 
CONDITIONS FOR THE MAJOR PLANETS OR USF THF SCHIJ- 
BAR T-S TUMP FF CONDI TONS FOR THE SOL AR SYSTEM A T 
FPOCH J. D„=?4?ooon. 5 . 

INI T=0.. .. .USF S TUMPFF— SCHURAR T INITIAL CONDITIONS. 

I N I T = 1 . . . . . US F R SUBMITS OWN INI TI AL CONDI TIOMS 
IMPLICIT REAL*8 (A-H,0-Z) 

REALMS A(20,3,17),R(?0*3, 12), HI 20,70) 

COMMON / I NOU T / IN , I OUT 

IN = 5 

I0UT=6 

READ (IN, 1 ) NS I Z E , I N I T 
FORMAT (215) 

CALL CONTRL ( A , R , D, MS I Z F , I N IT ) 

R E TURN 
END 

SUBROUTI NE COM TRL ( A , R , D , MS I ZE , I MI T ) 

IMPLICIT REAL*R (A-H,0-Z) 

*** HE I DFLBERG M-RODY PROGRAM 

BY J • SCHURAR T AMD P.S TUMPFF 

*** IRM S/360 FORTRAN 4 VERSION 

BY P. COMFLLA Z B. LOWRFY 

GODDARD SPACE FLIGHT CENTER 
GREENBELT, MARYLAND 20771 

ADAPTED FROM 

*** YALE FORTRAN 4 VERSION *** 

BY M, COOKE AND J. SCHUBART 


EXPLANATION of INPUT 

M ORDER OF THE INTEGRATION (M.LE.5).IM THE INITIAL 

ITERATION THIS CORRESPONDS TO 7*M DIFFERENCES, IN 
THE EXTRAPOLATION COMPONENTS TO 2*M+1 DIFFERENCES. 

I EXP 10**(-IEXP) LIMIT OF ACCURACY FOR THE INITIAL ITERATION 

N NUMBER OF BODIES N.LE.60 

KO =G*EM(1 ), WHERE G IS THE UNIVERSAL GRAVITATIONAL CONSTANT 

IN UNITS OF L***3/( (S 1 **)*^) 

E M ( I ) -MAS'S OF THE I TH BODY IN UNI TS OF M. 

MOTE: IF A BAR YCENTR I C SYSTEM IS NOT USED FOR INPUT, 

THF RELATIVE SYSTEM MUST ORIGINATE IN BODY *1. 

XP ( K , I ) ^SPATIAL COMPONENTS OF THE I TH BODY IN UNITS OF LENGTH, L • • , 
(IN THF APPROPRIATE CO-ORDINATE SYSTEM). 

XDOT( K , I ) 

=VELOCITY COMPONENTS OF THE I TH BODY IN UNITS OF L f,, /S», 
WHERE S IS THE UNIT OF TIMF. 

W =THF CONVERSION FACT OR NECESSARY TO EXPRESS KO IN UNITS 

OF <L**3/S**2 )**.5 

DIST = THE CONVERSION FACTOR NECESSARY TO EXPRESS XP IN 

UNITS OF L 

VFL =THF CONVERSION FACTOR NECESSARY TO EXPRESS XDOT IN 

UNITS OF LAS. 

H =THF INTEGRA HON STEP-LENGTH IN UNITS OF S 

T =THE STARTING EPOCH IN ARBITRARY TIMF UNITS. 

DEL TAT =THF INTEGRATION STEP-LENGTH IN THE SAME TIMF UNITS AS T. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 

c 

c 


( I PI lMCH d ) t 
I =1 t N ) 

I PUNCH ( J)=l 
= 2 
= 3 
=4 
= 5 
= 4 
=7 
= R 
=9 
= 10 
=11 
= 12 
= 13 

♦GE* 1^ 

IFG=0 
IFG.NF. 0 
K0=0 HR K0=2 
K0= 1 HR KD = 3 
IFD=Kn 

IF0=KD+4 


BARYt H*P . , PRINT 
RAR Yt 0. P * PUNCH 
RARY f O.P * tp.R I NT AND PUNCH 

R AR Y y S.P* vP R "I NT 

RARYvS.P.tPDMCH 

RAP Y* S* P • »PR INt AND PUNCH 
RELA f D.P. f PRINT 
REL A vU* P • -f PUNCH 
RELA,n.P, f PRINT AND PUNCH 
PFLA, S.P. * PR INT 
RFL A t S • P • t PUNCH 
RELAYS. P..t PR INT AND PUNCH 
NCI PRINT OUT 
NF W RUDY RFING ADDFD 

INPUT IN BARYCFN'TRIC SYSTFM. THI S SYSTEM MSFP fO». 

INPUT RELATIVE TO SYSTFM ORIGINATING I* 1 > 

nilTPUT IN RARYCFMTRIC SYSTFM, S.P. nR " •' P * * ^tJvfV Y 
niiTPUT TM RELATIVE SYSTFM, S.P. n R • O.P. RESPECTIVELY 
PRINT CPI— ORDINATES AND VELOCITIES FUR STFP Ml IMRFR 
AS GIVEN BY ILFM( I) , I=N7, N7-1 , , , , , , I 
PUNCH THIS OUTPUT AS WELL. 


INTEGRATION 


I LEM ( I ) , I = 1 , NT STEP NO. 


N7 


IZA 


AT WHICH TO PR I NT /PI INCH USING 
TO PRINT OUT CO-ORn. A NIT VEL 


FLFMMT. 
, DURING 


NO. OF TIMES 

Ei'rST^STEP AT WHICH TO PRINT OUT STFP NIIMRFR , T , 

• ;;^ H r^nZ^H STFP ciu^OROCKF TO PRINT/PUNCH 
COORDINATES OF BOOIFS ACCORDING TO VALUES OF I PUNCH 
I AST STFP AT WHICH TO PRINT THIS. 

WHEN STFP IZN IS COMPL. F TFO, 

NFW DATA ARE READ IN , FI THFR NEW ROOI F S OF MASS ZFRD OR 
Mfio I EYING CONDITIONS, THF PRES F N T C.A L C • IL A T 1 0 N S B E I w G 
r nM T T NU F n OR ENTIRELY NEW C AI.C.OL A TI ON-S ARF BEGUN, 

IE PRESENT CALCULATIOMs'aRE TO BE CONTINUED AFTER N P W BODY IS ADDED 
THE FOLLOWING MUST BE DONE (AT THE MINIMUM).--- 


IOFLTZ 


IZN 


SFT MASS, FM(J)=n, FOR EACH NFW RODY 
IF jmpuT COORDINATES and VELOCITIES RFLATIVE* 
for EACH NFW BODY e v , _ 

SFT I PUNCH ( J ) =X+13* WHERE 1. LE. X. LF, l . 
ACCORDING to DESCRIPTION ABOVE. 

INPUT V FI nr I TI ES AND coordinate, 
when 1 7 = 1 ZN GD TO ?:A00 new BODIES of ZERO mass 

r„"M C uI!z~ E r,S‘”Toon».T, 4 l .iz F «ss 

AND CO-ORDINATES FOR START OF NFW CASE. 
CONTINUATION CASE. CHFCK THAT ALL NFW BOOTFS OF 
7.FR0 HASS AMD INPUT THF IR CO—ORD I NA TFS . 

00 NOT COMPUTE ORRITAL ELEMENTS 
COM PI I TF ORRITAL ELFMFNTS 

=0 DO NOT COMPUTF ELFMFNTS FOR ROOY J 
=1 COMPUTF THEM 


A « 
R. 


C. 

I GO si 

I GO = 2 

I GO = 3 

inRR=n 

=1 


cnDMATfUHn f OOR DI-WAT FS OR MASSES INCOMPLETE) 

Sft r 5. «T Fn L Lnw»r, stpp W( >110) I 


20 



o-oooo 


1 1 CALL HR R I T PARAMFIER ' /7X, * J' ,5X, 1 IORRITI J) : 1 

2 'ORBIT I/O CONTROL P A R A m F T E R S 1 / ( I R , R X , I 2 ) ) 

1021 FORMA T ( * CALL ORIiCKE PARAMETERS : 1 /4X , 1 I Z A = ' ,16/ 

1 * I DEL TZ = ' , I6/4X, f IZM=» , I 6 / 7 X , ' J ' , 5 X , ' I PUWCHCJ ) * 

2 ' : DROCK E I/O CONTROL PARAMETER 1 / ( I R ,RX , I ? ) ) 

1022 FORMA T( 'OFLFMNT I/O PARAMETERS :' /4X, 1 I F0 = *, j ?/5X, 

1 ' KD = 1 1 1 2/ 1 CALL EL EMM T PARAMETERS: ' /5X, *N7=' , 14) 

10 23 FORMA T( ' 0 « , 5 X, 1 H= * , 01 ft . 8/6X , 1 T= 1 ,M ft. 8/ • OEL TA T= ' , 

1 016. 8 /5X, »K0=« , D?0.1 2 /ft X, 'W=' ,D16.R/3X, 

2 * 01 ST= * , 016.R/4X, ' VFL = 1 , 016, R/6X * ' N = * , 13/ 

3 4 X , 1 I E G = 1 , I 3 , 5 X , 2 0 A 4 / *0* , 2X , ' J V, T29 , * EM ( J ) • , 

4 T53, * X ( J) N T77, 1 X DO T ( J ) 1 / ( I 4 , 5X , 024. 1 ft , 5X , 024. 1 ft , 5 X , 

5 024.16/38X f 024,lft f 5X t 024,16/ 3BX, 024. 1 ft, 5X^024. lft ) ) 

10 24 FORMA T( '0* , 2X, ' J ' , T29, ' FM { J ) * , T53 , ' X ( J ) * , T77, ' XPUMKT ( J ) f / 

1 { I 4, 5X ,024,1 6,5X,D24.16,5X , D24 . 16/38X , 024. 16 , 5X , 

2 024.1 6/ 38X, 024. 1 ft , 5X , 024. 1 ft ) ) 

DIMENSION I PUNCH ( 50 } , I L E M ( 1 00 ) ,OIFMAX( 50) ,OIFJIIL ( SO) ,K000I 7.(50) ,01 
— SMI N ( 50 ) , 01 SjJUL ( 50 ) , KPD I S ( 50 ) 

DOUBLE PRECISION A ( MS I 7_ E , 3 , 1 2 ) , R ( NS I ZF , 3 , 12 ) , 0 { MSIZF, NS I 7 F ) 
SCHURART-STUMPFF INITIAL VALUES FOR F P OC H,J.O. =2430000. 5. ORDER 
OF INITIALIZATION AS FOLLOWS: SUN— MFRCt IR Y SYSTEM, VFNIIS, 

EARTH-MOON SYSTFM, MARS, JUPI TER , SA TURN, I JR AMI IS , NFPTUMF . PL I ITO. 
POSITIONAL COORDINATES RELATIVE TO CENTER OF SUN-MERCURY SYSTEM. 


IN A. 11, VS. 

REAL-8 XP ( 3 , 50 ) /3*O.DO f 

A -.5 11 3942059, -.4780976854,-. 1R030874R10, VENUS 

R -.261 498991 7, +.8696 237687, +.3771 652 157, E-M 

C -1.295477589, -.841 41 361 41, -.3513513446, MARS 

0 +3. 429472643, +3. 35 3869719, +1. 354948917 t JUPITER 

E +6. 641453441 00, +5. 971 56984400, +2. 18231501500, SATURN 

F +11 .263041 2500, +14. 6952588800 , +6 . 2796058 3300 , URANUS 

G -30. 1552293400, +1.65700086000, +1. 43785811000, NEPTUNE 

1 -21 .1238378000, +28. 4465 110 LOG, + 15. 3882665500, PLUTO 

H 1 2 3 * -1 . 2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 30- 2 / 

C VELOCITY COMPONENTS IN A.U./40 EPHEMFRIS DAYS. 

REAL *8 XD 0 T( 3, 50 1/3*0. DO , 

J +.56637681 8200, -.51 2087 158900, -.266497874500. VENUS 

K -.674657318300,-. 1700948 00800, -.07377431 92000, E-M 

L +.344004260500, -.369667484300, -.178937395200, MARS 

M -.222864773900, +.202 276882600, +.0922305 178000, JUPITER 

N -.1662? 88 59000, +.146271 2 35000, +.0676563647000, SATURN 

n -. 130 1308 16500, +.0758805 177900, +.0350896779200, URANUS 

P -.00961959898400, -.11 5065704000, -.046888 7522600, NEPTUNE 

2 — . 0707448500700, — . 0865592722000 , — . 00594685071 300, PLUTO 

0 1 23-— 1 .23456789101112130-2/ 

COMMON / TIME/ T, H, HO, HO, HI 0, HI 00, HOI 


COMMON /BLOCK 1/ CAI( 5,11 ), EM( 50 ) ,X(50,3) ,XPI1NKT( 50,3) , RFSCHL ( 50 , 3 ) 
1 , XN ARL A (50,3) 

COMMON / I M 0 F X / N , M t M 1 1 M 2 , M P M , M P 1 , I Z , I E G , I E n , I E X P , T P I J M C H 
REAL*4 I E G M S G { 2 0 , 2 ) 

DATA I EGMSG/ • I NP UT C 1 0-OR *,« 0 J NA * TF S * , 1 YS TF • , 1 M : R«, 

1 • ARYC * , ' ENTR * IC S10-* INP»,»UT C 1 , 'O-OP • , • 01 NA • , 

2 'TE S » , • YS TF 1 , *M: O 1 , • R I G I ' , * M A T F 1 , 1 S IN»,» ROD','Y *1 ' , 

3 8# ' */ 

COMMON /COMBLO/ GG( 11 ) ,00( 11 ) 

COMMON / I MO! IT/ I N, I OUT 

COMMON /ORR/ PI ,WK,0EG, IORR,IORRI T( 50) 

DATA SEVENS /-l . 234567891 01 1 1 2 1 30-2 / 

NAMELIST /INPUT/ N , FM , KO , W , XP, 01 ST , XOOT, VEL , T ,0FL TA J,H , M , 

1 I EG, I EXP, I ORR, I0RRIT,N7, I L FM , KO , I FO, I Z A , I DEL TZ , I Z.N , 
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2 I PUNCH, I GO 

REAL*R KO 
C 
C 

PI=DATAN2 ( 0 .00, -1 .00 ) 

M = 5 

IEG=1 

IFXP=14 

T=2430000.5 

deltat= 2 . no 

C UNIT OF INTEGRATION STEP 40 FPHFMFRIS HAYS 

h=.o25oo 
W=1 .00 

VFL=1 .no 
niST=i .no 

K Q = .473459612] 66R7n0 
0EG=180.n0/PI 

c 

c coefficients 
c 

00(1 ) =0.000 

D0( 2 )=0. 833333333333333330-1 
00 { 3) =0.833333333333333330-1 
00(4 ) =0.791 666666666666670-1 
00(5)=0. 750-1 

00( 6 ) =0.71345899470899471 0-1 
00( 7) =0.682043650793650790-1 
00 ( 8 ) =0.654957561728395060-1 
00(9) =0. 631 4043209876 54320-1 
00 ( 10 ) = . 61 072649861 71 2 3620-1 
00(11 ) =.592405641 233766230-1 
GG ( 1 ) =0.16666666666666667 
GG ( 2 ) =0.416666666666666670-1 
GG( 3) =0.22222 222222 2222220-1 
GG( 4) =0.145833333333333330-1 
GG( 5 ) =0.1061 50793650793650-1 
GG(6)=0. 8 22 585 9788 35978840-2 
GG( 7) =0.66479276895 9435630-2 
GG( 8 )=0. 55 37229938271 60490-2 
GG( 9) =0.47 180677475 81 6365 0-2 
GG( 10) =.409197067400192400-2 
GG( 11 )=. 359975497202679740-2 
I F ( IN I T.EO.O ) GO TO 2 
100 on 1 I = 1 , NS I Z F 
FM( I ) = S F V F N S 
oni j=i , 3 
XP( J, I )=SF\/FMS 

1 xnn T{ j , i )=sfvfms 

2 RFAO( 5, INPUT, FN0=1000 ) 

WRI TF( I OUT, 1020) I FXP, I ORR , ( J , I ORB I T( J ) , J=1,N) 

WRITE ( IOUT, 10? 1 ) I Z A , I OFL TZ , I ZN , ( J , I PI INCH ( J ) , J = 1 , N ) 

WRI TE ( IOUT, 1022) IF0,Kn,M7 

IF(M7.LE.O) GO TO 439 

WRI T E ( IOUT, 1017) ( ILFM( J), J=1,N?) 

N71 = l 

439 WRI TE ( I OUT, 10 23) H, T, DEL 1 A T, K 0 , w , n I ST, VFL , N , I EG , 

1 { IEGMSGi.!, IEG+1 ) , J = l, 20) , ( J, FM( J) . (XP( I J ) » X DO T ( T , J ), 

2 1=1,3) ,Jsl,N) 

MPM=M+M 

M1=MPM+1 
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M2-MPM+2 

MP1=M+1 

HQ*H*H 

HOI =* IDO^HO 

HD=l.DO/H 

H10=,in0*H 

H100= # 01D0*H 

on 4 1=1, N 

I F ( EMC I ) • EG • SEVENS ) GO TO 39 
I F ( I PUNCH ( I ) • L T. 14* AND* IGO, GT. 2 ) GO TO 4 

00 3 J = 1 f 3 

I F ( XP C J , I ) » F0« SEVENS* OR. XOO T( J , I ) • FQ* SEVENS ) GO TO 39 
X( I,J)=XP( J,I) *DI$T 

3 XPUNKTC I, J )=XDOT( J, I } *VEL 

4 CONTINUE 

I F { I GO • G T * 2 ) IG0=IG0-2 
WK=KO*W*W 

33 I F ( I Z A * L T# M ) IZA = M 

35 ICHFCK =0 

DO 405 1 = 2, N 

IF ( EM ( 1 ) *NE *WK ) E M ( I ) =WK * FM ( I ) / EM { 11 

400 I F ( I PUNCH ( I ) • L T* 1 4 ) 0,0 TO 405 

401 I E C EM ( I ) *NE • 0 * ) GO TO 406 

1 PUNCH ( I ) = I PUNC H ( I ) — 1 3 
ICHFCK = 1 

403 DO 404 K-1,.3 

X C I , K ) =X { 1 , K ) + X ( I , K ) 

404 XPUNKTC I , K ) =XPUNK T( 1 , K ) + XPUNK T ( I,K) 

405 CONTINUE 

I F ( I CHECK* EO. 1 ) IFG=0 
EM ( 1 ) =WK 
GO TO 40 

406 WRITE (6 ,1407) 

1407 FORMA T( 44H ADO I TI ON OF NEW RODY FORR I ODEN |N THIS CASF/ 

1 20H0 READY END 

RETURN 

39 WR I TE (6,1011 ) 

GO TO 100 

40 WRI TE C I OUT ,1024 ) { I , EMC I ) , ( X ( I , J ) , XPUNK T ( I » J ) , J = 1 , 3 ) » I =1 ,N) 

IF ( IORR.EO.l ) CALL ORRIT 
DO 440 J = 2,N 
OIFMAXCJ ) =0* 0D0 
DIFJULC J)=0.000 
KOODIZC J)=0 
D I S M I N ( J ) = 1 * 00+35 
DISJUL ( J ) =0*0 DO 

440 KPDIS ( J ) =0 
DO 441 J = 1,N 
DO 441 L = 1 * N 

441 DC J, L ) =1 *0035 
CALL KOEFFZ 
I Z = M 

CALL ANFITNC A,R t D,NSIZF) 

DO 45 1=1, M 
45 T= T+DELTAT 
GO TO 47 

46 CALL SCHRTTC A, R ,0, NS I ZF ) 

t=t+del tat 

IZ=IZ+1 

, DO 466 J=2,N 
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on 461 k = 1 * t 
ARR=OAR S ( R .{ ,K,M?) ) 

IF(ARR.LE.OIFMAX( J ) ) GO TO 461 

460 OIFMAXI J )=ARR 
niFJIIL ( J ) =T 

koooizc j)=k 

461 CONTINUE 
L F== J-l 

DO 463 L =1 , L F 

I F { n { j , l ) • G f . n i s m i n ( j ) ) an in 463 

46? OISMIN( J)sOU,L ) 

OISJUL ( J ) =T 
KPOI-St J')=100*L+-J 

463 CONTINUE 
IF(J.GF.N) an TO 47 

464 LA=J+1 

no 466 L = L A f M 

IF (0(L, J).GF.OISMIN(J) ) GO TO 466 

465 01 SMI N ( J ) =0 { L , J ) 

OISJUL ( J) = T 

K POI S ( J )=100*L+J 

466 CONTINUE 

47 IF (I Z — I Z A ) 51 *50*48 
4R IF ( IS-IOEL TZ ) 49, 50 , 50 

49 IS=IS+1 

an to 51 

50 IS = 1 

CALL ORUCKFI A,R,MSIZE) 

51 IF ( I Z • ECV. I ZN ) GO TO 53 

I F ( N7. L E. 0. OR. N7I • GT. N7 ) GO TO 46 
IFMZ.NE.ILEMCM71 ) ) GO TO 46 

53 CALL FLEMMT (B,NSIZE) 

WR I TF ( 6,1530 ) (OIFMAX ( J) ,01 FJUL ( J ) , KOOOI Z ( J ) ,01 SMIN{ J ) ,DISJULU) , 
1 KPni:.S< J),J=2,N) 

1530 F0RMAT(21H0 OIFMAX AMO OISMI M/ 

1 ( 1H , 1PE15.R, 1X,0PF10.1,1X, I1,5X, 1 PE 1 5. fi, 1 X , 0PF1 0. 1 , 1 X , I 4 ) ) 

N71=N71+1 
00 531 J = 2,N 
OIFMAX ( J) =0.000 

0 IF JUL ( J ) =0* 000 

Konoizi j)=o 

0 1 SMI N ( J ) = 1 . 0 0+3 5 
OlSJtlt ( J ) =0.000 

531 K PO I s ( J ) =0 

54 IFUZ.LT.IZM) GO TO 46 
GO TO (2,100,1000), I GO 

1000 RETURN 
FMO 

SUBROUTINE KOFFFZ 
IMPLICIT RFAL*R (A-H,0-Z> 

R E A L * R C A I { 5 ,1 1 ) , C ( 1 0 , 5 ) , R ( 1 1 , 1 1 ) , S 1 , S , T , T 1 

RF A L*B FAC (11) /l. '00, 1 . 0-0, 2. 00 ,6. 00,24.00, 1 * 20? , 7. 202,5.0403, 

1 4.03204, 3.6?RR05,3.62RR06/,SIG( 1 1 ) / 1 . 00 , -1 .00 , 1 . 00 ,-l .00, 

2 1 . 00 ,- 1 . 00 , 1 . 00 ,- 1 . 00 , 1 . 00 ,- 1 . 00 , 1 . 00 / 

COMMON /BLOCK 1/ CAI 

1 / I NOFX/ N , M, Ml , M? 

S=1.00 
on R L = 1 , M 

Li=?n 

L 2=L -1 
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SI =S**2 
on 7 k=i,li 
C(K,L)=0.00 
K 2 = K— 2 

I F ( K 2 ) 4,4,3 

3 C(K,L ) =C ( K , L ) +C ( K 2 , L 2 ) 

4 I F ( K + 1 -L 1 ) 5, ft, 7 

5 C ( K , L ) =C ( K , L ) -S 1 *C ( K, L 2 ) 
com 7 

ft C(K,L )=C(K,L )— SI 

7 CONTINUE 

8 s=s+i.no 
C L 1 =2*M+l 

L1=M1 
L2=L 1-1 
S-l.no - S 
nn in l=i,li 
R(Ll,L )=1. no 
nn 9 k = l , l 2 

K2-I.1-K 

9 R ( K 2 » L )=C{K2,M)+S*R(K2+1 ,L ) 

in s=s+i.no 
s=i .no 
no 13 l=i , m 
nn 12 l = i , L l 
L2=M2-L 
S1=S**2 
1 = 1 . no 

CAI( I, L )=o, no 
nn i i k = i,li 
T 1=R(K,L )*S1/T 
T=T+l.nO 

CA I ( I , L ) =C A I { I , L ) + T1 / T 

11 S1=S1*S 

12 CAK1,L)=CAI(I,|. )*Sir,( 1.21/ (FAC (l_)*FAC(L?> ) 

13 s=s+i.no 
RETURN 
ENn 

SUBROUTINE AMF I T M ( A , R , 0 , N S I 7 F ) 

AN FA MGS I TER A TI ON - STARTING ITERATION 
IMPLICIT RFAL*R <A-H,0-7) 

OIMFNSinN FMS(SO) 

nni lBLE PRFC ISION A ( NS I Z F , 3, 1 2 ) , R ( NS I z F , 3 , 1 2 ) , n ( n s T 7 F , MS I 7 F ) 
nnURLE PRECISION C A I , EM, X , XPUNK T, XNARL A » BFSCHI. » H * WO , FF , F I , EMMA , 

1 CAI0,Hl,H2,nuMMY 

COMMON /RLOCKl/ C A I ( 5 , 1 1 ) » FM ( 50 ) , X ( 50 . 3 ) , XPI INK T( 50,3) ,RPSCHI.( 50,3) 
1, XNARL A( 50, 3 ) 

COMMON /INDEX/ N,M,M2,M3,M4,M],IZ, I FG, IFO, IFXP. I PUNCH 
COMMON/ TI MF/ nilMMY , H, HO 

IPG = o , RARYCFNTRIC INPUT *=s=* ELSE INPUT RFLATJVF 

I F ( I EG.FQ.0 ) GO TO 105 

100 FMMA=0.0D0 

nn 101 j=i,n 

101 FMMA=FMMA+FM( J) 

nn 104 k=i . 3 
x ( l , K ) = 0 . ono 

XPUNKTd ,K 1=0.000 

nn 102 j= 2 ,n 
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X{ 1 ,K )=X( 1 ,K }+FM{ J )*x< J ,K ) 

1 0 3 XPIIMKT( 1,K ) -XPHMKTI 1,K )+FM(J )*XPIINKTU,K ) 

X ( 1 , K ) =— X ( .1 »K ) / F MMA 

XP1IWK T( 1 ,K } =-XPIIMKT( 1 K ) /FMMA 

nn 103 j = ?,m 

X( .J,K )=X< J,K )+X{ ] ,K ) 

103 XPIIWKK J f K )=XPUMKT< j,K )+XPUMKT( 1 ,K) 

104 CONTINUE 

10 5 COM T I NI IF 

CALL W W ( D t N S I Z F ) 

00 1 1=1,3 
on ] j=3,n 

A( J,I,13)=X< J,'J) 

1 A< J, I, Ml ) -R F SCHL ( J , I ) 

0FLTA=10.n0** ( -I FXP ) 
on 3 1 = P , w 

F MS m =OFL TA'M OARS ( RFSCHL { I , 1 ) )+OARS( RFSCHL! 1,7)1 +OARS ( RFSCHL { I • 3) 
*) ) 

IF ( fms ( i ) .of j. ,nn-?R ) on to 3 
3 fms ( 1 )=i.on-7fl 

3 continue 
mi 1=0 

r. M 7 = M 4 . M 1 

C M3=l +M2 

c M4=M+M 

ff = 0.000 
nn 31 mm=i t M 
3 ] FF = FF — 1 .000 

4 F T = F F 

1 T = 0 

nn 70 1 = 1 , m 7 

J 1 = I ~M1 
13 = 11 

IF(I1) 5,30,5 

5 .11= -II 

5 H1=FI*M 

T F ( mi ) . gt . 0 ) nn to r 
7 H3=n,5O0*Hl*Hl 
p nn 13 J = 7 , M 
on 13 K = 1 , 3 
TF(Mrj,F0,0) nn m 1.7 
9 X( J, K ) =0.00 
on 11 L = 1 , m 7 
L 1 =1 

IF f T3,OT.O) nn TO 11 

10 Li =M3-L1. 

11 X( vUK)=X( J,m+CAI (11,11 )*R( J,K,L) 

X ( J , K ) =HO*XU,K ) 

nn Tn 13 

17 x( J,K )=H2*A( J,K,«l ) 

13 X{J,K)=A(J,K,)7 1 +H 1 *XPI «NK T ( J , K J + X { ,J ,-K ) 

CALL wn(0,NSI?.F) 

I F ( mu . FO . o } nn to ift 
34 nn 17 j= 7 , n 

IFUT.GT.o) nn Tn ir 

15 OFL = OABS ( RFSCHL (J, 1 )-R( J, 1, I ) ) + DA RS ( RFSCHL ( J , 3 ) - R ( J , 2 , I ) )+OARS 
* ( R F SCHL { J , 3 ) — R { ,1 , 3 , I ) ) 

I F ( of 1 . . l E • F m s ( j ) ) nn Tn 17 

16 I T = 1 

17 CONTINUE 
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ift on 19 J = 2?N 
on 19 K = l?3 

19 A( J?K? I)=BESCHL( J,K) 

20 fi=fi+i.odo 
IFMt.EQ.O) on TP 23 

21 DP 22 J=2,N 
nn 22 k = i ? 3 
DO 22 L=1?M2 

22 Bl J,K f L)=A( JfK.t-L) 

MU = MIJ+1 

GO TP 4 

23 IF(MU.EQ.O) G'P TP 21 

24 WR I TE ( 6 ? 2 5 ) MU 

23 FORMA T(40H1 STAR TI NG ITERATION. NO. OF ITERATIONS = ,I3) 

1 2 = M — 1 
OP 30 J=2,N 
DP 30 K = 1 ? 3 
XNABL A { J?K )=0.0D0 
DO 27 L = 1 f'M2 
CA I D=CA I (M? L ) 

IF( I2.E0.0) GO TP 27 

26 CA I D=C A I D— C A I ( I 2 ? L ) 

27 XNABLA( J?K )=XNARLA< J , K )+C A ID*A ( J?K ,.L ) 

XNARLAfJ?K )=H*( XPUNKTC J?K ) +H* XNAR L A ( J « K ) ) 

DO 28 1 = 1? M4 

L1=M2-I 
DP 28 L = 1 » L 1 

28 A( J?K?L )=A( J?K?L+1 )-A( J?K,L ) 

DO 29 L=1?M2 

L1-M3-L 

29 B ( J ? K ? L ) =A ( J ? K ? L 1 ) 

30 B( J,K,M3)=0. DO 
RETURN 
END 

SUBROUTINE SCHRTT ( A ? B ? D, MS I ZE ) 

EXTR A PPL A T 1 0N$ SCHR ITT ?PRDN| ING‘=2* m +'1 - STEP PE INTEGRATION 
IMPLICIT REAL** ( A-H?P-Z ) 

DPI IB L E PRECISION A < M S I Z E ? 3 , 1 2 ) ? B ( N S 1 1 E ? 3 , 1 2 ) , 0 ( NS I ZE ? NSIZE) 

DOUBLE PREC ISIPM CA I,EM t X,XPUNKT,BESCHL?XNARLA ? DUMMY, H? HO? S fDDtGG 
COMMON /BLOCK l / C A I ( 5 ? 1 1 ) * EM < 50 > ?X ( 50 ? 3 ) ? XPUNK T ( 50 ? 3 ) ? BESCHL ( 50 ?3) 
1? XNABL A ( 50,3) 

COMMON /INDEX/ Ml ? M 2 

COMMON/ TIME/ DUMMY? H? HO 
COMMON /CONBLP/ GG ( 11 ) ? DD < 11 ■) 

DO 2 J=2,N 
DO 2 K = 1 ? 3 
S=0. ODO 
DO 1 L =2 ? Ml 

1 S=S+DD(L)*B(J?K,L+1) 

XNABL A ( J ?K ) =XNARL A ( J ?K )+H0* { B I J? K ? 1 ) +S ) 

2 X(J?K)=X(J?K) +X.N A R L A ( J ? K ) 

CALL WW(D?MSIZE) 

DO 4 J = 2 » N 
no 4 K = 1 ? 3 
A ( J ? K ? 1 ) =BE SCHL ( J ? K ) 

DO 3 L = 1 ? Ml 

3 A(J?K?L+1 )=A(J?-K?L) — R ( J ? K , L ) 

DP 4 1=1? M2 

\ B ( J ? K ? I )=A( J?K, I ) 


27 



RETURN 

FND 

SUR ROUT INF w w ( n , NS I Z E ) 

IMPLICIT REAL #8 (A-HtO-Z) 

REAL*8 DtNSIZF t MSIZE ) t 0X! 3) 

COMMON /BLOCK!/ C A I ( 5 , 1 1 ) t FM ( 50 ) , X ( 50 , 3 ) t XPI INK T! 50 , 3 ) , RESCHL ( 50 , 3 ) 
It X N A R L A ( 5 0 f 3 ) 

COMMON / I NOE X / N 
COMMON/ T IMF/ DUMMY tH 
no 2 k=i, 3 
x Cl t k >= o.no 

DO ]. J=2 1 N 

1 X C 1 ,K)=X(1 tK )+FM(.J-.)*Xf.J.,K ) 

2 X(1,K >=-X( 1,K }/FM{ 1 ) 

DO 14 J = 1 » N 

DO 3 K = 1 f 3 

3 BFSCHL ( J * K ) = 0. 00 

no 14 1=1 *N 
I F { I — J ) 4 1 1 4, 7 

4 I F ( F M < I ) )5tl4t5 

5 A = F M { I ) /D U t J ) 

DO 6 K = 1 9 3 

6 RE SC HI. { J r K ) =R E SCHL ( JtK )+A^( X(I f K )-X{ J f K) ) 

GO TO 14 

7 I F ( F M ( I ) ) 9 1 8 » 9 

8 IF(FM( J) )9, 14,9 

9 DO 10 K = 1 T 3 

10 DX(K )=X( I,K )-X! J ? K) 

A = OX ( 1 )*D.X( 1 )+DX( 2 )*DX{ ? >+DX( 3}*0X< 3 ) 

D! I, J )=DSORT( A) 

D( J » I ) =A*D( I ♦ J.) 

I F { J-l )llt 14, 11 

11 IE ( FM ( I) 112,14,12 

12 A = E M ( I )/Q( J, I ) 

DO 13 K = 1 , 3 

13 BFSCHL ! J , K ) -BFSCHL { J,K )+A*DX(K) 

14 CONTINUE 
RETURN 
END 

SURROUTINF FLFMNT (B,MSIZE) 

IMPLICIT REAL*R (A-HtO-Z) 

C 

C COORDINATES AND VELOCITY* IFD=KD PRINT ONLY , IED = KD+4 WITH PUNCH 

C 

999 FORMA T ( / /5 1H COORDINATES AND VELOCITIES FOR STEP NUMBER 

1 1 6/9H0 T = 1 P D 17*9) 

1003 FORMAT! 3 I A4, H , A4 , I 2 * A4 , 1 PD22 . 1 5 , A4/ ) , 

1 3 { 2 A 4, 11, AA, I2,A4tlPD22*15,A4/) ) 

1004 FORMAT ( 32H00 .01 #H#MAX (LAST USED DIFF.) = ,1PF10.2 ) 

1005 FORMA T ( 1 4H0 RODY MUMRER , I2,6H EM= ■- , 1 PD23. 1 5 ) 

1006 FORMA T ( 9H COORD. = , 1 P023 . 1 5 , 6X , 9HVEL DC I TY= , 1 PD?3 . 1 5 ) 

1007 FORMA T ( 9H COORD. = , 1 P3D16. 7, 1 1H VELOCITY* ,1P3D16.7) 

100R FORMAT! A4, 1PD22.15,A4) 

DOUBLE PRECISION R(MSIZF,3. 12) 

DOHRLF PRECISION CAI f FM, X , XPUNKTtBESCHL , XNARLA ,GG , Dn,V( 6) ,5, T, 

1 H, HO,HD t H10vH100 

real*r MU 

REAL*4 ALPHA! 7) 

COMMON /BLOCK!/ CAI! 5, 11 ) , EM ( 50 ) , X ( 50 , 3 ) vXPUNK T( 50 , 3 ) * RFSCHL ( 50 , 3 ) 
It XNABL A ( 50,3) 


28 



COMMON / I NO F X / N,M,M1 ,M2,MPM,MP1 ,IZ , IEG, IFO, IEXP, I PUNCH! 50) 

COMMON/ TIME/ T,H,HO,Hh,H10,H100 
common /CDMRLn/ GG< 1.1 ) *00! 1,1 ) 

COMMON /OPR/ PI, MU, RAD, I ORR , I OR R I T ! 50 ) 

data alpha/ 1 x p ( , }= xnni( T=, » / 

IF { IFn.GF,4)WRlTR( 7,1008) ALPHA! 6) , T, ALPHA (7) 

On 3 K-1,3 

XP! INK T { 1 , K ) =o.nno 

On 2 J = 2 , N 

s=o.ono 

nn ] L=l,Mi 

1 S=S+GG(L )*R.( J,K, L + l ) 

XPUNK T! J t K)=H.n-*XMARLAI J,K)+H*{R< J,K '1 )*fi .500-S > 

2 X P I INK T ( 1 , K ) =XPI INK T ( 1,K ) + FM ( J ) *XP! JMK T ( J ,K ) 

3 XPIINKT.C 1 »K )=-XPUNKT( 1 ,K ) / FM ( 1 ) 

WPITF! 6,999) IZ, T 

IF (MOn( 160,2 ) .FO V 0 ) GO TH 200 

nn mo j= 2 ,n 

WRITE! 6, 1005) J , E M ( J ) 

nn 5 0 K - 1 ,3 

X 3 = K + 3 

V ( K ) = X ( J , K j-XIUK) 

50 V ( K 3 ) = X PUNK T ( J , K ) -XPIJMKT! 1,K ) 

I F ( I F n . GE • 4 ) WR J TF ( 7 , 1 00 3 ) ( ALPHA ( 1 ) , K , A L P H A ( 2 ) , J , A L PHA < 3 ) , V ( K ) , 

1 ALPHA (7) ,K=1,3) , ( ALPHA (4) , ALPHA! 5) , K , A!. PHA ( 2 ) , J , AL PHA ! 3 ) , 

2 V ( K + 3 ) , A L P H A ( 7 ) , K = 1 , 3 ) 

I F ( K 0 . EO * 1 ) Gn Tn 75 

WRITE (A, 100 A) V( 1 ), V(4),V{2) , V! 5) ,V(3) ,V( A) 

Gn Tn 100 

75 WRITE (A, 1007) V 
100 CONTINUE 

Gn to 500 

200 nn 300 j = 1 , n 

WRITE (A, 1005) J,FM( V J) 

IFCKn.RO.l ) GO Tn 275 

IF! I FO.GE, 4 ) UP I TF! 7, 1003) ( ALPHA! 1 ) ,K, ALPHA! 2 ) , J , AL PHA ( 3 ) , X ( J , K ) , 

1 ALPHA! 7) ,K=1 ,3) , ! ALPHA! 4 ) , A L P H A ( 5 ) ,K, ALPHA! 2) , J , ALPHA (3) , 

2XPUNK.T! J,K),ALPHA(7),K = 1,3) 

WRITE! A, 100A) ( X { J , K ) , XPUNK T < J , K ) , K = 1 , 3 ) 

Gn TO 300 

275 WRITE! A, 1007) (X! J,K ) , K = 1 , 3 ) , ( XPUNK T ( J , K ) , K = 1 , 3 ) 

300 CUNTTNUF 
500 RM = 0» 

nn aoo j = 2 , m 
nn aoo k = 1,3 
BN-OAR S ( B ! J,K,M2) ) 

AOO BM-OMAXl { RM ? RN ) 

RM = RM*H100 

WRITE! 6,1004) RM 

IF f inRR.EO.l ) C AiLL ORBIT 

RF TURN 

FMn 

SURRnUTINE nRRIT 

C EQUATIONS FOUMn IN NASA nnCUMFN T NO. X-643-67-198 

C RY RnRERT RL AMCHARn 

IMPLICIT RFAL*R ( A-H,0-7. ) 

REAL*R MU,P(4),TMCL 

COMMON /BLOCK!/ C ( 1 05 ) , X ( 50 , 3 ) , XHOT ( 50 , 3 ) 

COMMON / I MOI IT/ IN, TOUT 
COMMON / I N 0 F X / N 
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COMMON /ORB/ PI T MI!,RAn T inRR,IORRlT( 50) 
WRITE! I OUT, 1 ) 

DO 1000 J = 1,N 

IF<IORRITUUEO,0)G0 TO 1000 
PER 100=0 • DO 
RA =0 , DO 

Pl= X ( J , 1 ) -X (If) ) 

P2 = X(J,2)-X(1,?) 

P3 = X! J,3)-X(l,3) 

P40=P1*P1+P?^P?+PB^P3 
P4=D$0RT! P40) 

vi = xnoT( J 1 1 ) — xnnT ( 1 , l ) 

V 2 = xnoT( J,2)-xnnTiU2) 

V3 = X no T ( J f 3 ) ~X DO T ( 1*3) 
V4G=V1*V1+V2*V2+V3*V3 
V4 = DSORT( V40) 

SORTMI l=D SQR T ( Ml t ) 

DZ “ P 1 # V 1+P2*V2+P3#V? 

D Z Ml ) = nZ/SORTMO 
A I N V = 2 • DO / P 4- V A 0 / Ml J 
AR A I NV-D ABS { A I NV ) 

IF (ARAINV.LF.KD-20) GO TO 1500 
50 SOA I NV-D SORT ( AR A I NV ) 

C Z = 1 • D 0— A I N V* P4 
SZ=DZMU*SQAINV 

c 

c semi-major axis-a 

c 

A- 1 * / A I MV 

S Z 2 = D Z M U*DZ Mf I * A I mi 

c 

C ECCENTRICITY -FCC 

C 

ECC-DSORT ( S 7.?+CZ*CZ ) 

C 

C ANGULAR MOMENT! IM/UM I*J MASS 

C 

HI = P2*V3-P3*V2 
H2= V1*P3-P1*V3 
H3= Pl*V.2 -P2*VI 
H4= DSORT ( HI *H.l +H2-*H2+H3*H3 ) 

T F R M A = { 1 . D 0 / P 4 - A I N V ) / F C C 
TERMB=DZ/ C ECC*MII> 

P( 1 )=TFRMA*P1~TFRMR*V1 
P-1.2)* TERMA*P2-TFRMB*V2 

P { 3 ) = TERMA*P3-TFRMR*\/3 
C 

C INCLINATION 

C 

COS I =H3/H4 

SI n i =dsortc i.no-cnsi4'Cnsi ) 
IFtcnsi.NE.o^no) on to 75 
I NCL =90 * DO 
GO TO 100 

75 INCL=DATAM{SINI/COSI ) 

C 

C RIGHT ASCENSION OF ASCENDING NODF 

C 

100 IF (SINI .FO.O.DOIGO TO 150 

HUM = H 4* SINI 


30 



SdMG=Hl/OUM 

C0MG=-H2/nUM 

OMEGA= 0ATAN2(SnMG,CnMG) 

c 

C ARGUMENT OF PER1GFF 

C 

ARGPER=DATAM2( P (3)/SINI f P1 1 >*CnMG+P(2) *.S'OMG ) 

C 

C MEAN MOTION 

C 

150 XMO T= S 0 R TMU* A R A I N V * SO A | N V 
C 

C MEAN AND ECCFNTRIC AMOMOLY 

C 

200 IFtA.LT.O.DOTGO TO 225 
EZ=DATAN2<SZ*CZ) 

Gfl TO 250 

225 EZ = DLOGUCZ + SZ)/ECC) 

250 XMA=EZ-SZ 

I F ( A. L T« 0* DO ) XMA=— XMA 
C 

C TRUE ANAMOLY 

c 

300 E2-*5D0*EZ 

EZ=FZ*RAD 

DUM = ( 1 * DO + ECC ) / ( 1 . DO-ECC ) 

IF<A»LT»0 .DO >G0 TO 325 
XP=OCOS( E2) 

Y=DSORT ( DUM ) *DS IN ( E2 ) 

GO TO 350 
325 EX1 -DEXP { E2 ) 

EX2 = 1 * D0/EX1 

Y-DSORTI -DUM ) * .500* ( EX 1 -EX 2 ) 

XP=.5*<EX1+EX2) 

350 TA = 0ATAN2(YfXP) *2.00 
C 

C TIME FROM PERIGEE 

C 

375 TPER=-XMA/XMOT 

C 

C BYPASS PERIOD AND APOFOCUS CALC IF HYPERBOLIC 

C 

IF (A.LT.O.DO) GO TO 900 
C 

C PERIOD 

C 

PER IQD= PI*A**1.5/S0RTMU 
C 

C APOFOCUS 

c 

RA=A*tl.D0+ECC) 

900 INCL= INCL*RAD 

OMEGA= OMFGA*RAD 
ARGPERs: ARGPER*R AD 
XMA= XMA*RAD 
EZ= EZ*RAD 
TA= TA*RAD 

WRI TE ( I OUTt 2 ) J, A,ECC, I NCL * OMEGA , TA , ARGPER , TPER 
1000 CONTINUE 

RETURN 
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1500 W R I T F ( I OUT , 3 ) .J,AINV 
GO T.n 1000 

1 FORMAT { » ORnnY NO.V, T13, 1 SFM.J -MAJOR AXIS* ,T33, 1 FCC FNTR I C I TY * , 

1153, 1 I W C L I M A T I OM * , T7 3 , 1 R • A . OF ASCENDING MOOR • , 

2 T98 , 1 TRUE AMAMOLY* ,T113,VARG. OF p FR IGF F 1 ) 

2 FORMA T( IR,2X,8n?0.12/T113, »TIMF FROM PFR IGFF » / 1 1 OX, 020. 12 ) 

3 FORMAT! •OSFMI-MAJOR AXIS OF ROOY N0.V,I3,* IS TOO L ARGF . l./A=«, 

1 020.12) 

end 

SOB ROUT INF 0R|1CKF{ A,R,NSIZF ) 

IMPLICIT RFAL*8 (A-H,0-Z) 

REAL *8 MU 
RFAL*4 ALPHA (7) 

DOUBLE PRECISION A { NS I ZF, 3, 1 2 ) , R ( NS I ZE, 3, 1 2 ) 

OOURLF PRECISION C A I , EM , X , XPUNK T,S{8),T,H ,HD , H10 , HI 00 , HOI ,HQ 

1 , og , on 

1 , RFSCHL , XNARL A 

COMMON /RL0CK1/ CAI ( 5, 11 ) , FM( 50) , X! 50, 3) , XPUNK T( SO, 3 ) 

1 , RFSCHL f 50,3) , XNARL A ( 50,3) 

COMMON /INDEX/ N, M, M2, Ml, MPM, MP) , I Z, I EG, I ED, I EXP, I PUNCH! 50) 

COMMON /TIME/ T,H,H0,Hn,H10,H100 r HQl 
COMMON / C 0 N R L O / GO ( 11), DO (11) 

COMMON /ORR/ PI ,MU,RAO f IORR, IORRI T< 50) 

1000 FORMA T( 1H0/14H0 STFP NIJMRER , I8,9H T ■= ,10017.9) 

1002 FORMAT (14H0 ROnY NUMBER ,12) 

1003 FORM A T ( 1 4H0 ROnY MUMRFR ,I2/9H COORO. = ,103028.15/ 

1 9 H V E L . = , 1 P 302 8.15) 

1004 FORM A T( 1 4H0 ROOY NUMBER , I2/9H COORO. = , IP 3020. 7/ 

1 9H V E L . - , 1 P3020 *7 ) 

1005 FORMAT! 2 (A4, I1,A4, 12, A4, 1 P01 5 . 7 , A 4 ) /A 4, I 1 , A4, I 2 , A4, 1 P01 5 . 7 , A4 , 2 A4 
1. II ,A^-,I2, A4,lP015.7,A4/2( 2A4,I1 ,A4, 12, A4,lPni5.7,A4) ) 

1008 FORMAT! 3 ( A4, II , A4, I 2 , A4, ] PD 2 2 . 1 5 , A4/ ) , 

1 3 ( 2 A 4 , 1 1 , A 4 , I 2 , A 4 , 1 P 0 2 2 • 1 5 , A 4 / ) ) 

1007 FORMA T ( 34H0 0. 1*H#H*MAX ( LAST USED OIFF.) = ,1PE10.2> 

1008 FORMAT! A4,lPn22. 15, A4) 

DATA ALPHA/ 1 XP(, )= XOOT f T=, '/ 

00 3 K = 1 , 3 
XPUNK T( 1 , K ) =0.00 

00 2 J=2,N 
0=0.00 

on 1 L = 1 , M2 

1 0=O + GG(L )*R( J,K,-L + 1 ) 

X pt ink T ( J , K ) =HO*XN A RL A ( J , K ) +H* ( R ( J , K , 1 ) *0 . 500-0 ) 

2 XPUNK T ( 1 , K ) =X PUNK T ( 1 , K ) +FM( J )*XPIJNK T ( J , K ) 

3 XPUNK TM ,K )= -XPUNK T ( 1 ,K )/EM { 1 ) 

IP = 0 

on 999 J = 1 , N 

1 = I PUNCH! J ) 

IF! I. F 0.13) CO TO 999 
IF! I.GF.7) GO TO 350 
IP=IP+1 

GO TO ! 50,100, 150,200,250,300), I 
50 IF! IP.FO.l ) WRITE {8, 1000) IZ,T 

WRI TF(6, 100 3 ) J, !X( j,K) ,K = 1,3) , ( XPUNKT! J,K ) ,K = 1 ,3) 

GO TO 999 

150 IF! IP.FO.l ) WRI TF { 8, 1000 ) IZ,T 

WRITE (6, 1003) J , ( X ! J , K ) , K = 1 , 3 ) , (XPUNKT! J,K) ,K = 1,3) 

100 IF! IP.EO. 1 ) WRI TE (7,1008) ALPHA ( 6 ) , T, ALPHA ( 2 ) 

WRITE (7, 1008) { ALPHA! 1) ,K, ALPHA! 2) ,J, ALPHA! 3) ,X ( J ,K } , ALPHA ( 7 ) , 

1 K = l,3), (ALPHA (4), ALPHA! 5), K, ALPHA (2-) ,J, ALPHA! 3), XPUNK T(J,K) , 



? ALPHA(7) ,K=1,3) 

00 in 99 9 

?nn iFdP.FO.i) W R I TF .( 5 , 1 OOO ) IZ,T 

MR I Tf ( 5,1004) J, ( X( j,K ).,K-1 , 3) , ( XPUMKTC J,K ),Ka.l,?) 

on rn 999 

*00 IFUP.FO.n WRTTF<5,]000) iz, t 

MR I TF < 6 , 1 00 4 ) .1 , ( X ( J , K ) , K = 1 , 3 ) ♦ ( X PI INK T ( J , K ) , K = 1 , R ) 

750 J F ( TP.F0.1 ) Ul R I T F { 7 • 1-008 ) A-L PHA ( 5 ) ,T,ALPHM?) 

WRITFjt 7, LOOS) ( ALPHA ( 1 ) ,K, ALPHA ( 7) , J,ALPHA( 3) ,X ( J,K ), ALPHA <7.) , 
] K = 1 * 3 ) , { A. L PH A { 4 ) , ALPHA ( 5 ) , K , A L PH A ( 7 ) , J , A l.-PHA ( 3 ) , XP1 INK T ( J , K ) , 

? ALPHA ( 7 ) , K =d ,3) 
on rn 999 

350 on 400 K = 1 ,3 

S ( K + ? ) = X PI INK T ( ,U K ) - XPI IMK 1(1»K) 

400 S { K ) = X { J , K ) -X ( 1 , K ) 

IP=IP+1 

Qn TO (999,999, 999,999,999, 399,450,500, 550,500,550,700 ) , I 

450 TF(IP.FO.l) HR I TF (5,10 00 ) 77, T 
MR I TF (5,1 003 ) J,S 
on Tn 999 

550 IFdP.FO.i) W R I TF ( 5 , \ 000 ) IZ,T 
M R I T F ( 5, 100 3 ) J, S 

500 IF ( IP.FO.i ) MR I TF ( 7 , 1-008 ) ALPHA (A) , T* ALPHA. ( ? ) 

MR I TF( 7, 1005) f ALPHA (1 ),K, ALPHA { 7 ) , J, ALPHA ( 3) ,S< K ) , ALPHA (7) , 

1 K = 1 , ^ , (ALPHA (4) , ALPHA ( 5) , K , ALPHA ( 7 ) , J , A L PHA ( 3 ) ,S(K+3), 

7 A L P H A ( 7 ) , K = 1 ,3) 

GO TO 999 

400 IFdP.FO.d MR I TF ( 5 , 1000 ) IZ,T 
MR I TF( 5,1004) J,S 
GO TO 999 

700 IFdP.FO.i ) MR I TF ( 5 , 1000 ) 17, T 

MRITF(5, 1004) .US 

550 IF ( IP.FO.I )MRITF(7,100P) ALPHA (5) , T, ALPHA ( 7 ) 

MR I TF { 7 , 1005 ) ( ALPHA ( 1 ) ,K, ALPHA ( 7 } , .1 , A L PH A ( 3 ) , 5 ( K ) ♦ ALPHA { 7 ) , 

1 K = ) ,3) , { ALPHA (4) t ALPHA ( 5 ) , K , ALPHA ( 7 ) , J , Al. PHA ( 3 ) ,S(« + 3) , 

? ALPHA ( 7 ) , K =1 , 3 ) 

99Q COMTTHHF 

IF ( IP.FO.O ) OH TO ]?00 
AMsO.no 
no noo j=?,n 
nn non k = i , 3 

linn BMsO.MAXl(8M f nA.8S<8( J»K f -Ml) ) ) 

RM-ssH-0 1 ❖PM 

MR I TF ( 5 , ] 007) AM 

1 ?O0 IF ( 1 0 R R .FO.n CALL OR B T T 
PFTIIRN 5 
FMD 
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T.THff. 

000 20 1 

100 fiiNPtn IEG=1, 

200 KD=3, I ED=3 , 

300 XP ( 1 , 1 ) =0, , XP (2,1 ) =0. 0, XP ( 3, 1 ) =0. , 

XDnT(i t i)=o.o,xnnT(2,i )= o.o,xdot{ 3,n=o., 

XP( 1,2) =-0.5 113942959, XP( 2, 2) =-.4780976854, XP{ 3, 2) =-0,1830874810, 

XDOT ( 1,2 ) =+0.5663768 182, XOOT( 2, 2) =-0.5 12087 1589, xnOTi 3, 2) =-.2664978745, 

XP( 1, 3)=-0. 2614989917, XP ( 2 , 3 ) =+0 . 86962 37687 , XP ( 3 , 3 ) =+0 , 8771652157, 

XODT( 1,3) =-0.67466731 8 3, XnflT( 2, 3) =-.170094800 8, XOnT( 3, 3) =-0.07 377431 92, 

XP( 1,4) =-1.295477589, XP{ 2, 4 ) =-0.841 41 361 41 ,XP( 3,4) =-0.35) 3513446, 

XOf)T( 1 , 4)=+0. 3440042605, XDOT( 2, 4)=-0. 3696674843, XnnT( 3, 4 )=-0. 1789373952 , 

XP ( 1, 5 )=+3. 429472643, XP( 2, 5) =3. 35 386971 9 , XP ( 3 , 5 ) =+ 1 . ^54948917, 

XnOK 1 ,5 ) =— 0 .2228647739, XDOT (2,5 ) =+0.2022768826 , XDOT ( 3,5 ) =+0.0922305 1 78 , 

XP( 1,6) =+6.641 453441, XP( 2,6 ) =+5. 971 569844, XP ( 3 , 6 ) =+? . 1 82 31 501 5 , 

XOnT( 1,6) =-0.1 6622 88590, X DOT (2, 6) =+0.1 4627 12 350, XnnT( 3, 6) =+0.0676568647, 

XP( 1,7) =+11. 263041 25, XP( 2, 7) =+14. 6952 5888, XP( 3, 7) =+6.279605 833, 

XOO T ( 1,7) =-0.1 301 308 165, XDOT ( 2, 7) =+0.07 588 051 779, XnDT< 3, 7) =+0.0350896779 2, 

XP( 1,8 ) =-30. 15522934, XP( 2, 8) =+1.657000860, XP( 3, 8) = + 1.437858 no, 

XDOT (1,8) =— 0 .0096195989R4,Xn0T( 2,8 ) =-0.1 150657040, XOOT( 3,8 )=-0. 04688 P 752 26, 

XP( 1,9 ) =-21. 12 383780, XP( 2,9)= +28 .4465 1 1 0 1 , XP ( 3, 9 ) =+15. 3 8826655 , 

2000 XDOT ( 1 ,9 )=-0. 07074485007, XOnT( 2,9) =-0.0865592722, X0OT( 3,9) =-0.0059468507 13, 

2100 FM( 1 )=1. ,EM(2 > = 2,451E-06, FM( 3 )=3.03591E-06,EM(4)=3.2 3?6F-07,PM( 5 )=9.547R6E-04, 

2200 EM ( 6 ) =2 .8558E— 04, E M ( 7 ) =4 . 3727E-05 , EM( 8 ) =5 . 1 775 9 F -05 , FM( 9 ) =2 .78F-06 , 

2300 K0=. 47345961 216687, T=0.,DFLTAT= 2 . , H=. 05 , M= 5 , M7= 5 , W= 1 . , 

ILEM ( 1 ) =10, ILEM( 2 ) =50, I LEM (3 ) =500 , I L EM{ 4 ) =1000, ILFM( 5 ) =1059, 

I0RB=1 , IORB I T ( 3 ) =1, I Z A=1059, I DEL T2 = l 0 59 , I ZN=5295 , N = 9 , IFXP=] 0, 

I DRB I T ( 1 ) =0 , I ORB 1 T ( 2 ) =1 , I ORBIT ( 4 ) = 1 , IORB I T ( 5 ) = 1 , I OR R J T ( 6 )=1 , I DRR I T( 7 ) = 1 , 

I ORB I T4 8 ) = 1 , IDRB I T( 9 ) =1 , IG0=1 , 

2800 IPUNCH=9*7, 41*13, £FND IHPOT DATA - ORBIT OF LOST CEK METEORITE - PART I 

0029 CARDS 

FIGURE 1 
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SIN Pin IGn = 3» T=?440590. , nFLTAT=-2. , H=-.0*5, H * 10, 

M7 = 5, ILF.M<l)=10 f I L F M { 2 ) = 2 0 , ILFM(3)=30, ILFM(4)=30, IL FM( 5 ) =100 , 

I Z A=5?95 « I DEL TZ=400 , !ZN=54750, 

EM ( 10 ) =0. ♦ XP ( 1 » 10 ) =+. 6 1 198278 , XP ( 2 ♦ 1 0 ) =-2 .0232007* XP ( 3* 10) =-.87761397, 
X0r»T( 1 ,10) = + .32 109876, XDHK 3, 10)= + . 11411804, XHOT < 3 * 10 )=+. 13002863 * 

I PUNCH ( 10 ) =20 * £EN0 XHHJT DATA. - ORBIT OF DOST CHI METEORITE, PART 2 

FIGURE 2 
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